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Abstract 

Gravitational wave production from bubble collisions was calculated in the early 
nineties using numerical simulations. In this paper, we present an alternative analytic 
estimate, relying on a different treatment of stochasticity. In our approach, we provide 
a model for the bubble velocity power spectrum, suitable for both detonations and 
deflagrations. From this, we derive the anisotropic stress and analytically solve the 
gravitational wave equation. We provide analytical formulae for the peak frequency 
and the shape of the spectrum which we compare with numerical estimates. In contrast 
to the previous analysis, we do not work in the envelope approximation. This paper 
focuses on a particular source of gravitational waves from phase transitions. In a 
companion article, we will add together the different sources of gravitational wave 
signals from phase transitions: bubble collisions, turbulence and magnetic fields and 
discuss the prospects for probing the electroweak phase transition at LISA. 



1 Introduction 



In the next decades, a new science will emerge from direct detection of gravitational radia- 
tion, that will open a qualitatively new way of probing the distant universe. Ground-based 
(LIGO [1] and VIRGO [2]) and space-based (LISA [3]) interferometers will reach the required 
sensitivity to detect many kinds of distant sources over a range of more than a million in 
frequency. Because gravitational waves (GW) penetrate all regions of time and space, with 
almost no attenuation, GW detectors can explore scales, epochs and new physical effects not 
accessible in any other way. 

Although the first GW detections will come from astrophysical processes, such as merging 
of black holes, another mission of GW astronomy will be to search for a stochastic background 
of GWs of primordial origin. An important mechanism for generating such a stochastic GW 
background is a relativistic first-order phase transition [4,5]. In a first-order phase transition, 
bubbles are nucleated, rapidly expand and collide. The free energy contained in the original 
vacuum is released and converted into thermal energy and kinetic energy of the bubble walls 
and the surrounding fluid. Most of the gravitational radiation comes from the final phase 
of the transition, from many-bubble collisions and the subsequent MHD turbulent cascades. 
The associated GW spectrum encodes information on the temperature of the universe at 
which the waves were emitted as well as on the strength of the transition. The characteristic 
frequency of the waves corresponds to the physics that produces them. For cosmological 
processes, this is close to the Hubble frequency, H ~ T^/Mpi. Once redshifted to today, this 
corresponds to 

/ ~ 1 mHz ^ . (1) 

100 GeV ^ ^ 

Remarkably, for transitions occuring near the electroweak epoch, / is in the frequency range 
covered by LISA (10~^ — 10~^ Hz). It is therefore very exciting that LISA could help 
probing the nature of the electroweak phase transition, and therefore provide information 
that is complementary to the Large Hadron Collider and the future International Linear 
Collider. 

Many types of new physics predict first-order phase transitions. Electroweak symmetry 
breaking in extensions of the Standard Model may be associated with a first-order phase 
transition (see for example Ref. [6]). Besides, the last decade has seen the emergence of 
the "landscape picture", following developments in String Theory. Strongly warped re- 
gions (throats) in higher- dimensional space-time are generic features in the string-theory 
landscape [7] and the phenomenological consequences are only starting to be explored (for 
instance through the prototype of Randall and Sundrum [8]). One interesting aspect is 
the cosmological evolution in these backgrounds. Thanks to holography and the AdS-CFT 
correspondence, a change in the 5-dimensional metric as the temperature decreases can be 
understood as a confining phase transition in the dual 4-dimensional gauge theory, and in 
models like [8], we typically expect first-order phase transitions at the TeV scale [9-13]. 
Finally, phenomena such as preheating at the end of inflation could share some common 
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features, as far as gravity wave emission is concerned, with the physics of first-order phase 
transitions [14-16]. 

The GW spectrum resulting from bubble collisions in first order phase transitions was 
computed in the early nineties [17-20]. It was realized ten years after the original calculation 
of [17-20] that turbulence in the plasma could be a significant source of GW in addition to 
bubble collisions [21,22]. Subsequently, the authors of [23] studied the GW signal due to 
a first order electroweak phase transition in the Minimal Supersymmetric Standard Model 
(MSSM) and its NMSSM extension. More recently, model-independent analysis for the 
detectability of GW with LISA [24,25], LIGO and BBO [25] were presented, relying on 
the formulae derived in [17-22]. The spectrum derived in Ref. [17-20] was estimated using 
numerical simulations, and no alternative calculation was performed afterwards. As argued 
above, we believe this is of high interest and it is time to revisit this question. In this paper, 
we present an analytical calculation of the stochastic GW background resulting from bubble 
collisions onlj0. Since bubble collisions take place in a thermal bath, and since we want to 
extend our treatment to deflagrations, we use the energy-momentum tensor of the relativistic 
fluid in the vicinity of the bubble wall as the GW source, rather than the energy-momentum 
of the scalar field. The result we find is comparable to that obtained by numerical simulations 
although the peak frequency is parametrically larger. 

A deterministic spherically symmetric expanding bubble does not produce gravitational 
radiation by itself. The reason is, that the transverse and traceless part of the energy 
momentum tensor for a radial deterministic distribution of the velocity field is identically 
zero (as we demonstrate in Appendix To produce a non-zero background of GW, one 
has to account for the fact that, towards the end of the phase transition, the collision of 
bubbles breaks spherical symmetry and leads to a non-zero tensor anisotropic stress. In the 
numerical simulations of Refs. [17, 19,20], this is accounted for by evaluating the transverse 
traceless component of an "incomplete" energy momentum tensor coming from the portion 
of bubble wall that remains uncoUided at a given time. This energy momentum tensor is not 
spherically symmetric and has a non-zero tensor anisotropic stress component. The total 
tensor anisotropic stress is obtained by summing all the contributions from single uncoUided 
bubble walls. Each simulation provides a given configuration of uncoUided bubble walls; 
bubble nucleation and coUision being random processes, the GW power spectrum is obtained 
by averaging the results of several simulations. This procedure is valid under the thin-wall 
"envelope" approximation, i.e. when the transition is strong and the bubble front evolves 
detonation. 

In the analytical evaluation which we present here, the situation is quite different. The 
GW production comes not only from the bubble wall, but from the entire fluid velocity 
profile in the vicinity of the phase discontinuity. If the non-zero fluid velocity shell contracts 

^Turbulent fluid motions triggered by bubble collisions together with magnetic fields are actually ad- 
ditional relevant sources for gravity waves from phase transitions. These effects have been reexamined 
recently [26-28] and since the subject is not closed, we will present revisited results from these contributions 
elsewhere [29]. 
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to a surface with vanishing thickness, no gravitational waves are produced. Therefore, we 
are not working in the thin wall approximation, and this is why we are able to apply our 
results also to the case of deflagrations. 

The paper is organized as follows. In Section [2] we review the general procedure for 
calculating the relic energy density stored in a stochastic background of gravitational waves. 
Section [3] describes our model of the GW source, the calculation of the bubble velocity 
power spectrum and the anisotropic stress power spectrum. In Section H] we define the time 
dependence of the phase transition parameters. In Section [5] the calculation of the GW 
spectrum, applicable both for detonations and deflagrations is presented. In Section El we 
make some comments on our analytical approach. In the last section we collect our final 
results and compare them with the existing formulae used in in the literature. Some technical 
aspects related to the calculation of the velocity power spectrum are collected in Appendix 
IBI Appendix [C] is a discussion on the behaviour of the small and large scale tails of the 
GW power spectrum. While the existing literature provides approximate expressions for the 
peak amplitude and peak frequency of the signal, there is no justification for the shape of the 
spectrum. Our analytical approach provides a rationale for it based on simple dimensional 
arguments. 



2 Gravitational wave power spectrum: general remarks 

Our goal is to estimate the gravitational wave energy density generated by bubbles during a 
first-order phase transition. This kind of cosmological source leads to a stochastic background 
of GW, which is isotropic, stationary, unpolarized and therefore characterized entirely by its 
frequency spectrum [30]. We consider a Friedmann universe with fiat spatial sections. The 
tensor metric perturbations are defined by 

ds'^ = a^[-drf + (% + 2hij)dx'dx^] . (2) 

The gravitational wave energy density is then given by 

^^^^^^ = MGa{^f ■ 

The over-dot denotes derivative with respect to conformal time and (...) denotes both time 
averaging over several periods of oscillation and ensemble average for a stochastic back- 
ground. The variables x and later also r denote comoving distances, rj and later r, C, denote 
comoving time. The density parameter is always scaled to today, ^xijf) = Px{v)/ Pcivo), 
where the index q indicates the present time. For relativistic species we have therefore 
^xiv) = ^x{Vo)/'^'^{v)'^ normalize a(?7o) = 1 and sometimes denote the present value 
of a density parameter simply by flxivo) = ^x] likewise, pc = pdvo)- 'H = d/a de- 
notes the conformal Hubble parameter. The radiation energy density today is taken to be 
aad(r?o)/i' = ^r.dh^ = 4.2 X 10-5 [31]. 
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We define the statistically homogeneous and isotropic gravitational wave energy density 
spectrum by 

{k,ik,v)h:^iq,v)) = 5ik - ci)\h\\k,v) , (4) 

where k is the comoving wave vector. The gravitational wave energy density, normalized to 
the critical energy density is: 

aa„(„) = ^^^l (5) 

where the factor (27r)~^ comes from the Fourier transform convention. We want to estimate 
the present day gravitational wave energy spectrum, in other words the gravitational wave 
energy density per logarithmic frequency interval. 



dink 



k^\h\'^{k,r]o) 



In an expanding radiation-dominated universe, hij(k, rf) is the solution of the wave equation 

2 • 

/i,j(k, r/) + -/iij(k, ri) + /c2/i,j(k, r/) = STiGa^{ri)Y[ij{\i, r]) . (7) 

Uij (k, 77) is the tensor part of the anisotropic stress, the transverse-traceless component of 
the energy momentum tensor that generates tensor perturbations hij of the metric: 

where Pij = 6ij — kikj is the transverse projector and T/m(k, 77) are the spatial components 
of the energy momentum tensor. As will be discussed in the next section, the anisotropic 
stress is a stochastic variable for the generation process under consideration. It accounts for 
the intrinsic randomness of bubble nucleation and collision. 

Our source of gravitational radiation is active for an interval of time corresponding to 
the duration of the phase transition, which is much shorter than one Hubble time [32,33]. 
We can therefore neglect the expansion of the universe while the source is still active, and 
rewrite Eq. ([7]) as 

h'-j{x) + hij{x) = ^^^^^* Uij{x) , (9) 

where x = ki], ' denotes derivative with respect to x and a* is the scale factor at the time 
of the phase transition. The dependence of hij(k,r)) on directions of the wave- vector enters 
only in the polarization of the wave and is irrelevant for our discussion. As will become clear 
at the end of this section, in Eq. (fT6!) . this is due to statistical homogeneity and isotropy of 
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the source. We assume that the source turns on at time r/in and turns off at time ?7fin- The 
solution of ([9]) is 

hij{x < Xfin) = — p-^ J dyg{x,y)Uij{y) , (10) 

where y = kr [r denotes conformal time) and Q = sin{x — y) is the Green function satisfying 
Q{x,x) = and Q'{x,x) = 1. Once the source is no longer active, we have to match the 
above solution with the solution of the free wave equation during radiation domination 

h'y(x) + h',^(x) + h„(x)^0 (11) 

h,(x > X,.) = 4. /'°'^-^'"' + . (12) 

The matching procedure gives the coefficients 

Bij = p * XfinJ dy sm{x{ir, - y)Uij{y) , 

Aij = — H T:^x&n dy cos{x{in - y)'nij{y) . (13) 

In order to simplify the equations, we neglect the first term in Aij which gives a subdominant 
contribution to the GW spectrum in the range of frequencies we are interested in. In fact, 
this term contributes in a sizable way only for modes larger than the horizon, 

XRn = krjfin < 1 , A; < l/?7fin ~ K , (14) 

where Ti^, denotes the conformal Hubble factor at the time of the phase transition, and is 
assumed to be constant from ?7in to rjan since the phase transition lasts for a time much 
shorter than one Hubble time. We will see that the GW spectrum grows very steeply at 
large scales (as k^) and peaks at a scale corresponding to the maximal size of the bubbles, 
which is typically much smaller than the horizon. Therefore, we are mainly interested in the 
sub-horizon part of the spectrum and in order to evaluate it we can safely neglect the term 
Bij/xfin- Using definition (jlj) and solution ([T2|) we finally find {z = kQ 

\h\Kx)\'' = ^ {{A,A^,) + {B.,B:^)) 



"fin 

k^ J 2x^ 



dy / dz cos{z — y)Il{k,y, z) (15) 



In the double integral above, we have combined the products of two Green's functions into 
the simpler term cos{z — y). Moreover, we have introduced the unequal time correlator of 
the tensor anisotropic stress in Fourier space, 

(H,,(k, r)H*.(q, C)) = 5(k - q)H(A:, kr, kQ . (16) 
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The delta function is due to the statistical homogeneity of the source, and because of statisti- 
cal isotropy the power spectrum of the anisotropic stress only depends on the wave number. 
Note that for the matching we have used the free wave propagation equation ffTTj) . which is 
valid in an expanding, radiation dominated universe with 0(77) oc 77. Hence, solution f[T^ for 
7] > 7]^, implicitly assumes that the number of relativistic degrees of freedom is constant. We 
come back to this issue in section [5l 

To summarize, in order to determine the spectrum of the gravitational radiation Eq. ([6]), 
we have to calculate the power spectrum of the anisotropic stress evaluated at different times. 
This requires computing the correlator of the energy momentum tensor. The next section 
is devoted to a calculation of Il{k,y, z) (Eq. (TH]). For this we need a model of the energy 
momentum tensor that sources the gravitational waves. 

3 Model of the GW source 

We now develop a model for the stochastic source of gravitational radiation. We are dealing 
with a cosmological first order phase transition taking place in a thermal bath [20,34]. The 
cosmic fluid of the initial metastable phase supercools until the nucleation of bubbles of the 
final phase can start. The initial high-temperature phase or false vacuum is typically but 
not necessarily the symmetric phase. However, in the remaining of the paper, we will use 
the term "symmetric" for the initial phase and "broken" for the final phase. The phase 
transition ends when the entire universe has been converted to the broken phase by bubble 
percolation. We are only interested in the last stages of bubble growth. Towards the end of 
the phase transition, the bubbles can be considered simply as spherical combustion fronts 
moving at constant velocity [35]. Any memory of the initial shape of the bubbles, driven by 
the scalar field dynamics, is lost and the problem can be reduced to a purely hydrodynam- 
ical description. The bubbles are modeled as spherically symmetric configurations of fluid 
velocity. The velocity field is a stochastic variable, following the intrinsic stochasticity of the 
nucleation process. 

3.1 Anisotropic stress power spectrum: general remarks 

Since we are interested only in the anisotropic stress, we start with the spatial, off-diagonal 
part of the energy momentum tensor of the cosmic fluid, quantifying the spatial components 
of the kinetic stress-energy tensor of a bubble configuration [20]: 

T«fe(x, r) = (p + p)— — . (17) 

1 — ti'^^x, t) 

V is the velocity of the fluid in the frame of the bubble center, and v = ||v||. We want 
to calculate the anisotropic stress power spectrum given in Eq. f|T6l) . In order to simplify 
the calculation, we neglect the spatial dependence of the fluid enthalpy density w = p + p 
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and of the gamma factor 7^ = 1/(1 — v"^). This assumption is necessary in order to be 
able to proceed analytically. It supposes that the only stochastic variables in the problem 
are the fluid velocity components Va{'x, r), and that the spatially dependent 7 factor can be 
approximated by 7(x) ~ (7) = 7. The consequences of this assumption cannot be quantified 
exactly. However, we know that (f^) varies smoothly from f|(rint/-R)^ to vj (see Eq. E71) 
where rint and R are defined in Eq. [26l A conservative choice is to always set (f ^) to its 
smallest value, and this is what we will do in Eqs. (155|1 and (1561) . Under these assumptions, 
we can write the Fourier transform, 

r„b(k,r) = , "^^^1 , [ £pva{k-p,r)v,{p,r) . (18) 



With this expression, the power spectrum of the energy momentum tensor involves the 
four-point function of the velocity distribution: 



(r,,(k,r)r;,(q,c)) = 



d^p j d'^h{va{k-p,T)vb{p,T)v^{(i-hX)vd(h,C)) ■ (19) 



There is in principle no reason why our stochastic velocity field should have a Gaussian 
distribution. However, we have to make some assumptions in order to calculate analytically 
the four-point function in the above expression. As one often does, we assume that Wick's 
theorem, which is strictly valid only for Gaussian random variables, gives a good enough 
approximation to the four-point function. It certainly gives a better estimate than, for 
example, the simple product of expectation values. Applying it we find 

X j d'p[Cac{p,TX)CM{\k-p\,rX) + CUp,rX)CU\^-p\,TX)] , (20) 

where 

Cacip, rX)= I d'rCacir, r, Oe^^'"- , ^.(r, r, C) = {va{^, r)t;,(x + r, ()) • (21) 



The correlation between point x and point y = x + r is a function of r only because of 
statistical homogeneity. 

In our approach, the ensemble average in Eq. ([3]) is now traced back into a correlator for 
the bubble velocities. This is where the stochasticity of the process is encoded. Since the 
velocity field is statistically homogeneous and isotropic, its power spectrum has the general 
form (see the next subsection 13.2. ip 

Cacip, r, C) = F{p, r, 06ac + G{p, r, Qpapc • (22) 
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The power spectrum of the tensor part of the anisotropic stress is calculated using the 
definition (fT6|) by applying the transverse traceless projector as in Eq. ([8]): 



(n,,(k, r)n*^.(q, 0) = n,6cd(T.,(k, r)T;,(q, C)) 

Vabcd = (^PiaPjb — (^) {^icPjd — ^PijPcd^ (q) (23) 

A somewhat lengthy calculation yields 

n(fc, r, C) = (i_^r|l)||"ff,,2(^)) / [^np)F{\k - p|) + 2(1 - /3^)F(p)G(|k - p|) 
+ 2(1 - \')G{p)F{\\, - p|) + (1 - - P')Gip)Gi\k - p|)] (24) 

where \ = k ■ p and P = k ■ k — p and we have suppressed the time variables r and C in 
and G. 

The problem is now reduced to the determination of the functions F and G which define 
the power spectrum of the fluid velocity via Eqs. fl21ll22p . For this, we need a model of the 
fluid velocity which we discuss next. 



3.2 Velocity profile of bubbles 

Since we are only interested in the last stage of the phase transition, we consider the hy- 
drodynamics of bubble growth at late times, when a steady state solution is reached. The 
bubble wall, in other words the combustion front where the phase transition is happening, is 
moving at constant velocity. In the hydrodynamical description of the combustion, the front 
is treated as a surface of discontinuity. Energy and momentum must be conserved across 
the front and all the entropy production is confined to it [34,35]. Elsewhere, the fluid is 
in a state of thermal equilibrium. The energy momentum tensor of the burnt (broken) and 
unburnt (symmetric) phases is simply that of two perfect fluids (see Eq. [T7Il . There are two 
kinds of solutions to this hydrodynamical problem, detonations and deflagrations. These are 
classified following the characteristics of the fluid flow in the rest frame of the combustion 
front [36,37]. 

In detonations, the incoming velocity of the symmetric phase fluid into the front is 
supersonic vi > Cg in the rest frame of the front. The outgoing velocity of the broken 
phase fluid out of the front can be supersonic V2 > Cs for weak detonations, or equal to the 
speed of sound for Jouguet detonations V2 = Cg. The case of strong detonations V2 < Cg is 
forbidden [36]. Although weak detonations are possible [38], in the following we concentrate 
for simplicity on the case of Jouguet detonations. This is the case analyzed in [20], for which 
the dynamics of the bubble growth is completely determined in terms of the phase transition 
strength. Since both fluid phases are relativistic, they have the sound speed Cg = l/Vs. 
In the rest frame of the bubble center, the velocity of the bubble front Vh is supersonic. 
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corresponding to = ^'i > c^: the symmetric phase fluid is therefore at rest, and the front 
is followed by a rarefaction wave in the broken phase fluid. The rarefaction wave brings the 
fluid motion to rest towards the center of the bubble. Near the detonation front, the broken 
phase fluid velocity V{ in the rest frame of the center of the bubble is simply given by the 
Lorentz transformation 

Vi = 25 

1 - V1V2 

The velocity proflle of the broken phase fluid for a Jouguet detonation has been studied in 
detail in Refs. [20,34,36] and is shown schematically in the flrst panel of Fig.[TJ As customary, 
we show the velocity proflle as a function of the parameter r/t. Here t = is the time of 
bubble nucleation and r denotes the distance form the bubble centeiH. We remind that this 
situation corresponds to the steady state solution at late times, long after the nucleation 
time. The velocity of the broken phase fluid goes to zero in the interior of the bubble at a 
distance from the center corresponding to Cgt [36]. 

For deflagrations, on the other hand, the incoming velocity of the symmetric phase fluid 
into the front is subsonic vi < Cg. The outgoing velocity of the broken phase fluid out of 
the front can be subsonic V2 < Cg for weak deflagrations, or equal to the speed of sound 
for Jouguet deflagrations V2 = Cg. Strong deflagrations are again impossible [36]. In the 
rest frame of the bubble center, the velocity of the bubble front corresponds in this case to 

= ^2- The front moves at subsonic velocity = V2 < Cg and is therefore preceded in 
the symmetric phase by a shock wave. Inside the combustion bubble (broken phase) and 
outside the shock wave (symmetric phase) the fluid is at rest; in between, the fluid moves 
outwards. In the case of planar deflagrations, it does so at constant velocity given by V{ 
(Eq. [35,36]. The qualitative features of the velocity proflle for a planar deflagration are 
shown in the second panel of Fig. [T]E|. 

For our analytic calculation, we want to simplify the real velocity proflle both for deto- 
nations and deflagrations. We assume a velocity proflle which grows linearly within a shell 
near the bubble wall, as shown in the last panel of Fig. [H We have normalized the veloc- 
ity proflle at the outer boundary (bubble wall or shock front) to the correct value Vf for 
detonations and deflagrations. This is because the biggest contribution to the GW energy 
density comes from the highest velocity region. Therefore, our approximated proflle does 
reproduce the most relevant feature as far as GW generation is concerned. The boundaries 
of the shell, deflned as Viat and Vout in Fig. [T] are left as free parameters, in order to allow for 
an approximated description of both detonations and deflagrations. 



^Throughout this section, for simphcity wc use a generic time variable t; we switch back to comoving 
time in paragraph 13.2.21 

■^In Appendix A of Ref. [20] it is shown that, in the case of spherical deflagrations, the velocity profile 
actually decreases between V2 and Wshock- We do not account for this behaviour here, since in any case we 
are forced to introduce an approximate form for the velocity profile. 
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Figure 1: This figure shows the quahtative profile of the velocity of the broken phase fluid 
in the frame of the bubble center, for detonations (top panel), planar deflagrations (middle 
panel) and the approximation given in Eq. (l26l) (bottom panel). The horizontal axis shows 
r/t where t denotes the time after bubble nucleation [t = 0) and r is the distance from the 
bubble center. 



The simplified profile is 

( +\ - j (^f/-^) - ^o)a for Tint = fint^ < |x - Xq] < -R = fout^ , fr,(^^ 

Va[^,t) - I Q Otherwise. ^^^^ 

Here xq is the position of the bubble center. The radii of the shell's inner and outer bound- 
aries are respectively r^nt = Vmtt and R = Voutt, where t is much later than the nucleation 
time t = 0. In the case of Jouguet detonations, the inner boundary is in the broken phase and 
corresponds to t^int = Cg. For deflagrations, it corresponds to the bubble wall, Wmt = f 2 = v^. 
The radius of the outer boundary is R = and for detonations it is associated to the 
bubble wall velocity fout = "^i = "fb, while for deflagrations to the shock front Vout = "^^shock- 
To summarize: 

Detonations Deflagrations 

f int = Cs Vint = f 2 = f b 

Vout = Vi = Vb . Vout = fshock • 
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The numerical values of Vi^t ? ^out and V{ will be crucial in determining the amplitude of the 
GW signal and will be discussed in Sections 15.21 and I5.3[ A schematic drawing of the bubble 
(or shock front) is given in Fig. [2l 



R=v„,„t 




Figure 2: A schematic drawing of the non-zero velocity region, corresponding to the bubble 
(for detonations) or to the shock front (for deflagrations). 



3.2.1 Velocity power spectrum 

Given the velocity profile, we can proceed to calculate the velocity power spectrum. We start 
by evaluating the two-point correlation function at equal time for fixed positions x and y 
defined in Eq. f l^T]) . The position of the bubble center xq defined in Eq. fl2Bl) is the stochastic 
variable. Therefore, in the region of non-zero velocity we have: 



where we remind that V{ is the maximal value of the fluid velocity in the rest frame of the 
bubble center. For the velocity correlation function not to be zero, x and y must be separated 
by a distance |x — y| < 2R and they have to be in the same bubble (or shock wave, in the 
case of deflagrations). Moreover, they have to be in the shell where the fluid velocity fl2^ 
is not zero. These conditions are satisfied provided that the center xq of the bubble they 
belong to is in a volume Vt given by the intersection of two shells centered in x and y, which 
have inner radius Tint and outer radius R (see Fig. [3]). Therefore, the correlation function is 
given by the mean over all the possible center positions xq within this intersection volume, 
multiphed by the probability a{t) that there actually is a bubble center in this region: 



{vi{:>^,t)vj{y,t)) 




((x-xo)i(y- xo)j) 



(27) 



vll_ 




'xo(x - xo)i(y - xo)j . 



(28) 
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As customary in cosmology, here we use the ergodic assumption: ensemble averages are 
equivalent to space averages. The probability of having a bubble center in the intersection 
region is simply given by 

a(t)=0(t)| (29) 

where (j){t) is the fraction of volume occupied by bubbles at time t, and Vc is the volume of 
the region where xq can be, in order for x or y to be in the same bubble: the total volume 
of the two overlapping spheres in Fig. [31 Setting r = x — y we find 



3 V 2 8 



(30) 



The tensorial structure of the two-point correlation function of a statistically homoge- 
neous and isotropic field is known: the correlation function can only depend on the distance 
between x and y. We choose an orthonormal basis with x — y || 62. The off-diagonal 
components of the integral in Eq. (128!) are zero by symmetry, and therefore we find 



/jj(r, R, Tint) 



Iij{r, R, Tint) 



/ xo(x-xo)i(y-xo) 

(l)it) ^ y Iij{r, R, Tint) 

111 Sij + {I22 - hi) Tifj 



(31) 

(32) 
(33) 



(since In = 133). The functions In and I22 have to be calculated by performing the necessary 
integration in the four volume regions Vi shown in Fig. 131 The details of the calculations are 
given in Appendix B. 

The velocity power spectrum is then obtained by Fourier transforming the two point 
correlation function (l32ll with respect to the variable r. Remembering the definitions fl21|22l) 
one finds the general expressions (JF denotes the Fourier transform) 



{v,{k,t)vUq,t)) = 6{k-q)C,,{k,t) = 6{k-q)[F{k,t)6,,+G{k,t)kk,] 



F{k,t) 
G{k,t) 



„,2 r 



^'1 



1 d ^ 
JF 

k dk 



'22 



'11 



l± ( I22 - In 
kdk \ r^Vr 



dk^ 



I22 - 1 
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r^Vr 



where we remind that (j){t) is defined below Eq. (j29l) . 

We define the new dimensionless variable K = kR and the fraction 



(34) 
(35) 

(36) 



Vint/ Vout — l^int/ R ■ 



(37) 



Note that 1 — s = {R — r\at)/R is the relative thickness of the shell, and will contribute to 
the amplitude of the GW signal. In our approach, the GW signal will vanish in the limit of 
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Figure 3: This figure shows how the intersection volume Vi changes as a function of the 
separation between x and y, r = |x — y|, where x and y are located at the centers of the 
shells. The upper left, upper right, lower left and lower right plots respectively correspond 
to < r < R - Tint, -R - Tint < r < 2rint, 2rint < r < R + r^^t and R + ri^t < r < 2R. 
Therefore, this figure does not depict bubble collision (in our approach we do not actually 
collide bubbles). The shaded volume does not represent the volume of intersection between 
two different shells, but it accounts for all possible positions of the center of the bubble to 
which two given points x and y belong. 

vanishing thickness. We perform the Fourier transform and obtain the following expression 
for the velocity power spectrum: 



{v,{k,t)v*{q,t)) = 5{k - q) A7i(P{t)vfR{tY[A{K)5,j + B{K)k,k,] 



(38) 




K < 4.5 
K > 4.5 , 



(39) 



B{K) 




exp (-(0.7- 2.5)72) (^)^ if K < 0.7 

exp (-(A' - 2.5)2/2) ' if 0.7 < K < 4.5 (40) 

exp (-(4.5- 2.5)72) (^)^ if K > 4.5 . 
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The last two equations are good fits to the real functions A{K), B{K) which are shown in 
Fig. m A{K) behaves as white noise for K ^1 while B{K) oc K^. On small scales (large 
values of K) both A{K) and B{K) decay like K''^. In the region where A{K) and B{K) 
are maximal, our approximations overestimate the real functions for values of s > 0.65 (we 
will comment on this in Section [52]) • In the following analysis we will consider the maximal 
value of s = 0.74: for this value the approximations overestimate by about 16%. 

The behaviour of A{K) and B{K) can be predicted from general considerations. First 
of all, the correlation function fl28l) vanishes for r > 2R : the characteristic scale of the 
correlation function is the diameter of the bubbles. We therefore expect this scale to show 
up in the power spectrum at wave-number k ^ 27r/2i?. This is indeed what happens, since 
the functions A{K) and B{K) change their behaviour at approximately K ^ 2.5 ^ tt. 
Moreover, since the correlation function is a function with compact support, its Fourier 
transform, the power spectrum, must be analytic in k. We remind that an analytic function 
can be developed in power series around any point of its domain. The term A{K)6ij of 
Eq. (1391) is analytic for — s> 0, if and only if A{K) oc and n is an even integer n > 0. 
This justifies the white noise behaviour observed at large scales for A{K). On the other 
hand, analyticity of the term B{K)kikj for K ^ is satisfied if and only if B{K) oc K" 
with n an even integer and n > 2. This justifies why B{K) increases as at large scales. 

3.2.2 Unequal time correlation function 

Up to now we have evaluated the velocity correlation function at equal times. We note 
however, from Eq. (|2T|1 . that we actually need the correlation function evaluated at different 
(comoving) times r, The velocity in point x at time r can be correlated with the velocity 
in point y at time (. Consider, for example, t < (. In this case, the unequal time correlation 
function is not zero if the velocity shell in the bubble includes x at time r, and grows to 
include y at time (. According to the approach outlined above, evaluating the correlation 
function at different times means performing the volume integral of Eq. fl28|) within regions 
Vi given by the intersection of spheres of different radii (cf. Fig. [3]). This integral is too 
complicated to be done analytically: therefore, within our analytical approach, we first try 
to simply approximate the unequal time correlation function with the one at equal time 
calculated in the previous subsection. This is a reasonable approximation, provided that the 
region of non-zero velocity at r overlaps with the region of non-zero velocity at (. If this is 
not the case, we simply set the unequal time correlation function to zero. 

Reintroducing rj^^ as the time of nucleation we find that, in the limiting case, the inner 
boundary of the non-zero velocity shell in the bubble at time ( equals the outer boundary 
at time t ii ( = {t — rjinj/s + rji^, where we remind that s = rint/R- Symmetrizing among 
the two times, we have then 



~ (i;,(x, r)vj{y, r))0(C - r)e((r - r/i„)/s + - C) + 
+ {Vi{^, CMy, C))0(r - C)e((C - Vin)/s + Vm - t) 



(41) 
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Figure 4: Velocity power spectrum. The top left panel shows the function A{K) determining 
the diagonal part of the velocity power spectrum and the fit given in Eq. fl39|) for different 
values of s = fint/'^out = fint/R- The solid lines from top to bottom (red, yellow and pink) 
are the correct functions and the dashed lines (green, blue and cyan) are the fits for s = 0, 
0.6 and 0.75 respectively. The top right panel shows the function B{K) and the fit given 
in Eq. ( HOl) . again for the same values of s. The approximations overestimate for s > 0.74 
by about 16%. The lower panel shows again A{K) and B{K) and the fits of Eqs. fl39| HOl) 
for s = 0.6. The flatter curve is A{K) and the dashed line is its fit given in Eq. (l39l) . We 
note the white noise behaviour of A{K) for small values of K. The more peaked solid line is 
B{K) and the dashed line the fit given in Eq. fHUl) . At small K, B[K) grows like while 
A{K) is constant. At large K both functions decay like 
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where G(r) is the Heaviside function. We choose arbitrarily to set the time appearing in the 
equal time correlator corresponding to the smaller of the two times. 

In Section [5l we derive the gravitational wave spectra obtained using this approximate 
form of the unequal time correlation function and discuss its shortcomings. We will even- 
tually propose another method, which consists in giving an approximate form directly for 
the unequal time anisotropic stress power spectrum, rather than for the velocity correlation 
function. As we will see, proceeding in this way we can have better control over the positivity 
of the power spectrum, and obtain more reliable results. 

3.3 Anisotropic stress power spectrum: calculation 

We now have everything we need to evaluate the anisotropic stress power spectrum of our 
source, using Eq. (l24l) . 

The unequal time correlation function Eq. fl4ip . together with the equal time velocity 
power spectrum given in Eq. (1551) . lead to 

n(^, = (1 _ .^'^l) jl'ff - Rir? G(c - r)e (i(r - v..) + - c) 

X J £P [4A(P, t)A{\K - P|, r) + 2A(P, t)B{\K - P|, r)(l - A^) 

+ 2B{P,t)A{\K - P|,r)(l - P') + B{P,t) B{\K - P|,r)(l - - X')] 

+ symmetric t , (42) 

with P = pR(t), K = kR(T), X = k ■ p, and (3 = k ■ k — p. Using the functions A[K) and 
B{K) given in Eqs. ( l39l l40l) . we can perform the above integral. A good approximation to 
the integral, after factorizing out the s-dependence as (1 — s^)^, is given by 

1+ f-)' 

I(K) = 0.0412 ^ ^ . (43) 

i+(f)^+(fr 

This fit, together with the exact integral, is shown in Fig. [51 The function is flat on large 
scales and changes slope at ~ 3. A white noise behavior at large scales is expected, since 
the anisotropic stress power spectrum is the convolution of the velocity power spectrum: this 
simply means that the anisotropic stress is not correlated at distances larger than the source 
correlation scale 2R. In Appendix [Cl we derive this in details. As one sees in Fig. |5l Eq. (143!) 
is a very good approximation to the numerical integral. On small scales the convolution 
decays like A{K) and B{K), hence like This can be understood as follows: when 

K ^ -ft'max; wherc -ft'max dcuotcs the wave number at which A{K) and B{K) peak, the main 
contribution to the convolution integral comes from the region |P — K| ~ i^^max ^ K. The 
value of a typical term in this region is about A{K)A{K^s,^), and the phase space volume 
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is -R'max- hence we expect T{K) ~ K^^^A[K^a_^) A[K) for K ^ A'max and analogous for the 
contributions containing B{K). 

Before we are ready to insert the expression P2l) for n(A;, r, in Eq. f|T5|) to evaluate the 
gravitational radiation power spectrum, we need to determine the time dependence of R{'r]) 
and 0(77). 




K K 



Figure 5: These two figures show the exact integral in Eq. (H2l) in solid (black), and the fit 
I{K) given in Eq. ( 143|) in dashed (red). In the right panel we clearly see the white noise 
behavior for small values of K and the behavior for large values of K, as expected (see 
discussion above and in Appendix [C]) . 



4 Time dependence of the phase transition parameters 

We now investigate the actual time dependence of some parameters introduced previously, 
such as the fraction of volume occupied by bubbles at time 77, 0(?7), which we need in Eq. ( l29l) . 
or the bubble radiu^. In this section we closely follow Ref. [33] in the modeling of the first 
order phase transition. 

The rate of bubble nucleation of the broken phase bubbles is defined as r(?7) = M.^a^ e~^^'^\ 
where Ai is the energy scale of the phase transition and S^r]) the tunneling action. We 
Taylor expand the action at first order around a fixed time rifi^: the time at which the 
transition ends. Defining /3 = —dS / dri\rf^^, one can rewrite the nucleation rate as r(?7) = 

"'in this section we specify what we actuahy take for the bubble radius, and we will denote it by R{rj). We 
remind that, in the case of deflagrations, this does not coincide with i?, which is the position of the shock 
front, but with rint- 
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r(?7fin) exp (/?(?7 — r/fin))- Note that we define (3 = a^P in terms of comoving time, differently 
from tlie usual convention. The probability that a given point remains in the false vacuum 
at time t] is given by 

p{n) = e-'^"^^ (44) 

where I{ri) is the fraction of volume occupied by broken phase bubbles at time rj without 
considering bubble overlap [33,39,40]. Assuming that the universe remains static for the 
entire duration of the phase transition, and assuming a constant velocity for the bubble 
expansion V]^, I{ri) is simply given by 

m = - drV{r)vl{r^-rfc^STi^V{r^) (45) 



The quantity 0(?7) in Eq. (129|) is given by 0(?7) = 1 —piji). The times ?7in and ?7fin are defined 
such that p(?7in) — 1 and p(?7fin) — 0. Following [19,33], we choose a number M ^ 1 and 
define r/gn as r(?7fin) = /^^M/Svrwj^, so that piji^-n) = exp (— M) ^ 0. In the same way, we 
choose a number m <^ 1 such that [3{r|f^^ — rji^) = In (M/m) and p(77m) = exp (— m) ~ 1. 
This gives the duration of the phase transition 

van - Vin = In — . (46) 
m 

In order to evaluate the mean bubble radius, we consider the number of bubbles which 
have a given radius 6 at time rj. Calling rjs the nucleation time of a bubble with radius 6 at 
time rj, one has: 



dN 
~d5 



^iV5)p{vs) ^^^^ 



The shape of this distribution is shown in Fig. 5 of Ref. [33]. For each rj, it has a maximum 
at the value R{ri) = ^ ln/(?7): this value defines the mean radius of the bubbles at time 77. 
Calling f] the time at which p{f]) = 1/e, so that /(f/) = 1, we set 



-JO for r/in < r/ < r/ , 

R[V)-<y ^ln(/(r7)) foTn<V<mn. 



The condition I(f]) = 1 defines f] = rjaa — 1^ InM. 

As already explained in Sec. 13.21 bubbles can be treated as combustion fronts moving at 
constant velocity only at times much later than nucleation time. Therefore, we identify ?7in = 
f] in the evaluation of the emitted gravitational radiation. This leads to ?7fin = ?7in + In M. 
We further decide to neglect the logarithms and simply set 

r/fi,-r/i„ ~ (48) 
R{vi) ~ fb(^?-^in) • (49) 
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If we do not neglect the logarithms, we have to replace the final bubble size Vb(3~^ by 
ln(M/m). Neglecting the logs, we identify the duration of the phase transition with 
and the radius of the bubbles at time r] with its mean value Rit]). We do not account 
for the possibility of having bubbles of different sizes at a given time. 



5 Evaluation of the gravitational wave spectrum 

We now want to evaluate the gravitational wave energy density per logarithmic frequency 
interval ([6]) , which is given in terms of the GW power spectrum (fTS!) . According to Eq. (fTSjl , 
the latter evolves as oc ?7~^; however, this behavior is strictly valid only in a radiation 
dominated universe with a constant number of relativistic degrees of freedom {c.f. Sec. [2]) 
and thus lacks generality. To estimate the GW power spectrum at 77 > 77*, we therefore 
evaluate it at the end of the phase transition: \h'\'^{k,ri^,), for which no assumptions have 
been made concerning the number of relativistic degrees of freedom. Then, we simply use 
its radiation-like evolution: 



dQ{k,ri) dQ{k,r]^, 



dink dink 
From Eqs. ([6]) and f|T5l) . reminding that x = krj we have 

dn{k,r],) ^ k^\h'Wk,7],) 
dink ~ 2{2TTfGp^al 



(50) 



(51) 



\h'\^{k,x^) = - , * / dy dzcos{z -y)U{k,y,z) (52) 



2 V A;2 



where in the last equality we have set xgn — x^,. In the above equation, we further have to 
substitute expression fH2|) for the anisotropic stress source. As already discussed in Sec. El 
since the source is active for an amount of time much shorter than one Hubble time, we 
neglect the expansion of the universe while gravitational waves are produced. Therefore, in 
Eq. fH2]) we set the enthalpy density to a constant, denoted by = w{t) ~ u^(C)- Moreover, 
we eliminate the time dependence of the 7 = \/l — v"^ factors by substituting v{t), v{() with 
the constant fiuid velocity sv^, corresponding to the fiuid velocity of the inner boundary of 
the shell (c/. Eq. (1261) and Fig. [2]), remembering that s = ri^t/R- We explain the reasons for 
this choice below. We remind that the double integration in Eq. (l52l) is in time, with the 
notation y = kr, z = k(. With X(A') given in Eq. fll5]) . we finally obtain using Eq. flSUj) 



^iii *^ ^i 



X 



'-^(r)i?^(r)e(C - r)e ( -(r - rn^) + r/jn - C ) T{kR{T)) + symmetric y ^ z 



(53) 
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Let us first investigate the pre-factor in the above expression. The enthalpy density is 



^* = 7TP;ad P:ad= - (54) 

3 \9*J at 

where p*^^ denotes the radiation energy density in the universe, go = 3.36 and g^: denote 
the effective number of relativistic degrees of freedom today and at the time of the phase 
transition respectively. We also define a dimensionless parameter estimating the amount of 
kinetic energy present in the source, with respect to the radiation energy density. From the 
definition (fT7|) of the energy-momentum tensor: 

_ 4 ^ {sv{)^ , , 

^kin - s^^-'^i-^sv.y ' ^^^^ 

fi;ad p:ad si-{sv,y ^ ' 

We have chosen to define the above parameter in terms of the fluid velocity at the inner 
boundary of the velocity shell svf, which always satisfies SV{ < Cg both for detonations and 
deflagrations. This ensures that fikin/f^*ad < 1 to remain consistent with an isotropic FRW 
universe. 

Summarizing, we can rearrange the pre-factor in Eq. fl53l) in terms of the above defined 
parameters, of the conformal Hubble factor Ti.^ = a^H^ and of ^lra.d, as: 

The gravitational wave energy density is therefore proportional to the square of the kinetic 
energy density of the source, as one would expect. 

In the double integral of Eq. fl53l) . we define the new integration variable u = kR{y) = 
Vont{y — Xin), and the new dimensionless parameter 

Z=^. (58) 

P 

Following the discussion in Sec. HI the function 0(r) becomes 

0(r) = 1 - exp (- exp (1 + /3(r - r]Rn))) • (59) 
We finally obtain for the gravitational wave energy density spectrum today 

dn{k,r]o)h^ ^ 3 fgoY' ,2/^kin^v^*^'(l-^''' 

' 'rad*^ 



dink 2iT3\gJ V^*ad/ \f3j 



+ f du (1 - e-'''^^"'f u'l{u) sin f ) |> • (60) 

JsZ V "^out , 
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We recover the known result, that the amplitude of the gravitational wave spectrum is 
proportional to the square of the ratio between the duration of the phase transition and the 
Hubble time, Ti^,/ (3 [4,41]. If the parameter s = 1, this means that the fluid motions vanish 
everywhere, and therefore no gravitational radiation is produced. The divergence in s — > 
is only apparent, due to our deflnition of the kinetic energy density parameter in Eq. fl56l) . 

Expression ( l60l) gives in all generality the power spectrum of gravitational radiation 
produced by spherically symmetric, stochastic conflgurations of fluid motions (such as broken 
phase bubbles in a phase transition), characterized by a velocity distribution which increases 
linearly in a shell of inner radius rj^t = fint^ and outer radius R = fout^- The integral 
determines the shape of the spectrum, and depends on the values of s and of the velocity fout- 
These parameters should be chosen according to the physical situation under consideration. 

We expect the large scale part of the GW spectrum to increase as k^. As explained 
in details in Appendix C, this is a simple consequence of the fact that the source has a 
flnite correlation scale (corresponding in our case to the length-scale 2R). It is therefore 
a generic behaviour for causal sources. Conversely, the small scale part of the spectrum 
depends on the details of the source correlation function. This part of the spectrum is in 
principle affected by our choice of a linear growth for the velocity proflle in the non-zero 
velocity shell (cf. Sec. 13.21) . If we could account for the correct velocity proflle, the power 
law dependence of the GW power spectrum at large k might be different from what we flnd. 
On the other hand, the frequency at which the GW power spectrum peaks can again be 
predicted by general considerations. In fact, the velocity power spectrum given in Eq. fl38l) 
has a characteristic wave-number corresponding to the bubble diameter, /c ^ vr/i? ^ 2.5/ R. 
Also the anisotropic stress power spectrum changes slope at about k ^ ir/R ~ 3/R. For 
the gravitational radiation power spectrum we therefore expect a peak at approximatively 
the same wave-number k ~ 7r/i?(r7fin) given by the mean bubble size at the end of the phase 
transition. 

5.1 Unequal time approximations 

The expression (1601) has been derived using the approximation Eq. (1411) for the unequal 
time correlation function of the velocity fleld. Under this approximation, the anisotropic 
stress power spectrum at different times takes the form given in Eq. (H2|) . which has been 
substituted in Eq. (l52l) in order to evaluate the GW power spectrum. According to the 
deflnition of the GW spectrum |/i'(x)p, the time integral appearing in Eq. (1521) should give 
a positive result. Therefore, the unequal time anisotropic stress power spectrum should be, 
by deflnition, a positive kernel, such that 

/ dy dzgix,y)g*ix,z)U{k,y,z)>0 , (61) 

Jxin JXin 

where g{x, y) generically denotes the Green function of the wave equation iQ. Even though 
the approximated form for the unequal time correlation function for the velocity fleld Eq. (14T!) 
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seems reasonable from a physical point of view, it does not lead to a positive kernel for 
n(fc, I/, z) as it should. Expression (l60l) is not always positive and therefore it is unacceptable. 
To avoid this problem, we now define approximations directly for the unequal time correlator 
Il{y, z) which are positive by construction. We first discuss two cases which are physically 
less motivated, but are very useful for comparison: the completely incoherent and completely 
coherent approximations. Finally, we discuss a third and better motivated approximation for 
n(?/, z) that still gives a positive result. We believe that this last approximation is the one 
that best recovers the true result. Apart from being strictly positive, it has the advantage 
to be very close to the physically motivated approximation for the unequal time correlators 
of the velocity field. As we shall see, the peak position of the GW spectra is very similar in 
all cases, while the peak amplitude varies by up to one order of magnitude. 

We consider first a totally incoherent source, that is uncorrelated for unequal times t ^ C,. 
We take the Ansatz (cf. Eq. |i2D 



(n,,(k,r)n^(q,C)) = 5(k-q)n(A;,r,r)^^^l^ (62) 



n(A;,r,r) = * .,., [^Mr)vn'R{rf{l - s'niKir)) . (63) 

We multiply the 6- function by the characteristic time 1/(3 in order to maintain the correct 
dimensions. Under this assumption, the oscillating term in Eq. (1521) . cos(x — y), integrated 
with the delta function, becomes simply 1 and we recover a positive result. Using the same 
notations as in Eq. (pUj) . the GW spectrum from totally incoherent bubbles is (see Fig. E]) 



dink 



incoh. 



s(i+"/z)n2 3 



X -/ du{l-e-''^''yu'Iiu). (64) 







The opposite situation is given by a totally coherent source, that is correlated for every 
time r and (. In this case, the Ansatz for the anisotropic stress power spectrum is 



(n,,(k, r)n^(q, 0) = 5(k - q) v!i(^v/n(^ 

V^Kk^) = -^^4vr0(r)^f2/2(r)3/2(l - s')VTiKW) (65) 

We rewrite the oscillating term cos(x — y) using the duplication formula. This allows us to 
split the double integral into a sum of two positive terms. The GW spectrum is in this case 
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(see Fig. [6]) 
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^rad 
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^'out 
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sm 



fout 



(66) 



The positive part of the resuh based on the approximated unequal time correlation 
function for the velocity field, is in between these two extreme cases: the Heaviside functions 
introduce correlation only among times r and C, which are sufficiently close to each other. 
We therefore expect also the correct result for the GW spectrum to be in between the 
limiting cases described above. Moreover, for sufficiently large scales the details of the time 
correlation do not matter, and we expect that the GW spectra derived using these different 
approximations become comparable at these scales. As we will see below, this is indeed the 
case (c./. Fig. [6]). 

Finally we introduce top-hat unequal time correlations directly in the anisotropic stress 
power spectrum, rather than in the velocity correlation function. This approximation is 
also an intermediate case between the completely coherent and incoherent ones, and it is 
constructed to always give a positive result for the GW power spectrum as we will discuss. 
We make the following symmetric Ansatz: 



(n,,(k,r)n*,(q,c)) 



n(fc,r) 



(5(k - q) [n(A;, T)Q{kC - kT)Q{x^ - {kC - kr)) 
U{k, CMkr - kCMxc - (kr - kQ)] , 



[SVi 



\2\2 



(67) 
(68) 



Under this assumption, we set the correlation to zero for modes with a time separation 
larger than Xc/k, where positive, dimensionless parameter of order unity (to be 

specified later). Therefore, the anisotropic stresses at different times are correlated if the 
time separation is less than about one wavelength. Physically, this just means that longer 
wavelengths are correlated over a longer time. This corrects the lack of correlation that 
resulted from our very first approximation. The GW power spectrum becomes: 
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) ^ 'rad 
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du(l — e 
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dv cos 
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(69) 
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which is apparently positive for Xc < vr/2. Resolving the Heaviside functions, we have to 
consider two separate cases. For Z < XcVout we find the spectrum (see Fig. |6]) 



dink 



mix 



^ 'rad I'' 



2Wout , /-I _e{i+i'/z)N2 3—' ^ . fZ — U 



/ du{l- e-^'-'^'Y u'I{u) sin ) , ^ < x.^out , (70) 

^ Jo V "^out 



and for Z > Xct'out we find 



dhik 



s4 



sin(xe) / rf«(l-e-'^""^"V«'X(«) 



+ r du{l-e--'''^"'"'fu'l{u)sm(^^)\, Z > x,v,^, . (71) 

JZ-XcVout V "^OUt / J 

The GW spectrum remains positive for < a^c < vr since in the range Z < Xc^out, sin {^f-f) > 
0, and in the range Z > XcVont, both sin(xc) and sin {^f-f) are positive. In the limit Xc — > 
the result tends to zero. A reasonable value is 7r/2 < Xc < vr. In Fig. [7] we show how the 
power spectrum depends on the value of Xc- We now present the results obtained for the 
different approximations discussed above, in the cases of detonations and deflagrations. 



5.2 GW from detonations 



As explained in Sec. 13. 2[ in the case of detonations R{t) = Voutt is the outer radius of the 
bubbles. Therefore, Vout = v\y = vi. We restrict ourselves to the case of Jouguet detonations, 
so that Viat = Cs = 1/ \/3- The value of the bubble wall velocity in Jouguet detonations is 
given in Refs. [20, 34] in terms of the ratio 

a = Pvac/P*ad (72) 

where p*^^ denotes the radiation energy density in the universe: 

Cs + \/(y^ + 2a/3 
= YV-a • ^^'^ 

The outer maximal fluid velocity is given by the Lorentz transformation (see Eq. [25!) 

f b - \pi{\/'io? + 2a - a) 

■'^'f = :; = : , (74) 

2 + 3a - V3a2 + 2a 
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where for the last equahty we have substituted fl73l) and = 1/ v^. Moreover, the parameter 
■s = t'int/'Wout takes the form 

s = ^ = . (75) 

For detonation^, it is therefore sufficient to specify a, and all the quantities necessary to 
evaluate the integrals in Eqs. fl60|64ll66ll70|7ip are determined. Knowing in addition the 
duration of the phase transition 7i*//5 fully determines the gravitational wave signal. We 
choose a broad range of values for a: a = 0.1, 1/2, 1, 10. a > 0.1 corresponds to s < 0.74. As 
discussed in section 13.2. H for s = 0.74 our approximations for the velocity power spectrum 
overestimate the true result by about 16%. 

Let us first discuss the shape of the spectrum. We have four different approximations for 
the unequal time correlator: the top-hat unequal time correlator for the velocity field in real 
space, Eq. ( !60|) . the incoherent case Eq. (IMI) . the coherent case Eq. (1661) . and the top-hat 
unequal time correlator for the anisotropic stress in Fourier space, Eqs. fl70ll7ip . They give 
rise to comparable spectra. In Fig. E] we show the result of the integrals as a function of 
Z = kvont/P- We have fixed the values a = 1/2 and Xc = O.Ovr. The spectra increase like at 
large scales, and have comparable amplitudes. This is expected, since the details of the time 
correlations should not matter at sufficiently large scales and the source is uncorrelated in real 
space (c/. Appendix C). Also, the positions of the peak approximatively coincide in the four 
cases, and correspond to ^peak — 4.6, so /cpcak — 4.6/3/wout = 4.6/i?(?7fin) — 7r/i?(?7fin)- This 
is comparable to the peak of the anisotropic stress power spectrum ~ 3/ R{t) , due to the 
fact that GW production accounts for the full evolution of the bubbles, and R{t) < -R(?7fin) 
for all times r < r]fi^. The peak of the GW spectrum is independent of time and corresponds 
to the mean bubble radius at the final stages of the phase transition, close to i?(?7fin). The 
amplitude at the peak is roughly the same for the two top-hat approximations and for the 
coherent case, with the value of the integral at the peak of the order of 0.08 ± 0.02. Within 
the precision of our analytical evaluation this difference is negligible. On the other hand, 
the incoherent case overestimates the amplitude by a factor 5. The small scale part of 
the spectrum is also different in the incoherent case with respect to the others: it decays 
slower, as as opposed to Z-f^ with (3 = 2 ± 0.2. The reason for this behaviour is 

apparent by looking at the integral in Eq. (l64ll : once Z has overcome the value at which the 
integrand peaks, the contribution from the integral becomes small and the function decays 
nearly like Z~^. Conversely, in the cases of Eqs. fl66|70ll71|] . at sufficiently large scales the 
spectrum scales nearly like Z"^. We explain the small scale decay on the basis of dimensional 
arguments in Appendix C. The two top-hat approximations are in very good agreement up 

■'To be in the detonation regime means that there is essentially no friction. In other words, interactions 
between the bubble wall (the Higgs field) and the particles in the thermal bath can be neglected. This is a 
strong assumption which can work if the released latent heat is indeed very large. Model-independent studies 
(such as Ref. [25]) have used Eq. ([75)1 for the bubble wall velocity. However, when considering a particular 
model, one should compute friction effects and therefore derive a more realistic value for the bubble wall 
velocity, following for instance the procedure of Ref. [42]. This requires to compute the bubble wall profile. 
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Figure 6: Integrals determining the GW spectra as a function of Z = kvout//3 for different 
approximations, in the detonation case with a = 1/2. The top (green), bottom (blue, oscil- 
lating) and middle (red) curves are respectively the incoherent (Eq. [61]), coherent (Eq. [66!) 
and top-hat in wave number (Eqs. fl70|7ip with Xc = O.dir) approximations. The middle 
line with the spikes (black) is the absolute value of the approximation with top-hat in the 
velocity correlation (Eq. [SU]). The spikes represent the passages through zero of this un- 
physical spectrum. All the spectra are comparable at large scales and have a similar peak 
frequency. The incoherent approximation overestimates the peak amplitude by nearly an 
order of magnitude. 



to Z ~ 10. For larger Z, the spectrum obtained for the top-hat in the velocity correlation 
becomes negative. This is due to the contribution from the first integral in Eq. ([60]) . which is 
negative after sZ > 8. Its absolute value is between the totally coherent approximation and 
the top-hat in wave-number approximation. We consider the top-hat in wave-number ansatz 
which is given in Eqs. fl70ll7ip to be the best approximation for the unequal time correlators. 

In the top left panel of Fig. [7] we show that the integral of Eqs. fl70ll7ip is only weakly 
dependent on a. The parameter a plays a greater role in determining the overall amplitude 
of the GW signal (which will be discussed later on). On the other hand, the shape of the 
spectrum depends significantly on the choice of Xc once Xc approaches vr. We remind that 
Xc defines the time interval Xc/ k beyond which the correlator of the anistropic stress tensor 
vanishes. In the top right panel of Fig. [7| we fix a = 1/2 and vary Xc- Values of Z around 
the peak correspond to the region Z > Xcfout, where the result is dominated, for a wide 
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Figure 7: The top left panel shows the weak «-dependence of our result Eqs. ( ITOfTTI) : the 
different (black, red, blue and green) curves are respectively for for a = 0.1, 1/2, 1, 10. The 
top right panel shows the dependence on Xc, for fixed a = 1/2. We remind that Xc defines 
the time interval Xc/k beyond which the correlator of the anistropic stress tensor vanishes. 
The lines from top to bottom (black, red green, blue and magenta) plot the integral in 
Eqs. ( 17011711) evaluated at Xc = 7r/2, 0.97r, O.lvr, 0.997r, 0.9997r respectively. The lower panel 
shows the GW spectrum of Eqs.f l70|l7T]) and its approximation Eq. (!76!) for a = 1/2 and 
Xc = 0.9n. 



range of values of Xc, by the first integral in Eq. flTTl) . The pre-factor of this integral decays 
like at high enough values of Z. This is roughly the decay we see in Fig. [7] for low Xc- 
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When Xc becomes very close to vr, the factor sin(xc) multiplying the first integral in Eq. (17T1) 
becomes so small that the contribution from the second integral in (17T!) takes over. For 
high values of Z, the second integral decays much faster than However, because of the 
integration limits, it is more and more suppressed as Z increases: therefore, at some given 
Z value, the first integral takes over again and the spectrum decays roughly like Z^"^. More 
precisely we find a Z"^ *^ decay at large Z. This behavior is reached for higher and higher 
values of Z as Xc approaches vr, as is shown in Fig. [71 For our final results, we choose the 
value Xc = O.dir. This may underestimate the signal somewhat at high frequencies compared 
to values x^ ~ vr/2, but this is a reasonable conservative choice. Values Xc still closer to tt, 
would be unjustifiably fine-tuned. 

The following approximate gravitational wave power spectrum is fairly general. We chose 
a = 1/2 but as shown in Fig. [71 the shape of the spectrum is almost insensitive to the value 
of a; the a dependence is essentially only in the prefactor (^kin/^rad)^ implicitly in 
(1 - s'Y/s^ : 
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where int|deto(^) = ; — ri — - — -jj , = 3.8 . (76) 
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This approximation for int(Z') is shown in the lower panel of Fig. [TJ and we remind that 
^kin/^rad defined as a function of s and Vf in Eq. f[56l) . 

Let us now investigate the overall amplitude of the signal at the peak frequency. 
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where the factor 0.084 accounts for the contribution of the integral at the peak frequency 
k ~ 4.6/?/uout (c/. low panel of Fig. [7|). In order for the bubbles to percolate and convert the 
entire universe to the broken phase, the phase transition must occur much faster than one 
Hubble time. We study the values of p/H^ ~ C(10), O(IOO), C(IOOO). Since the velocity 
of the bubbles = t'out is fully determined by the parameter a, so is the fluid velocity V{ 
and, in turn, the factor s = Cg/vout and the mean kinetic energy parameter ^tin/^vady 
Eqs. f [73|74|75p . We plot s and ^tiri/^rad ^ function of a in Fig. [8] Specifying a and /3 
fully determines the amplitude. Substituting in ([77]), and setting a = ^/3a^ + 2a, we find: 



dink 



1.7-10' 



{a - a)^{-3a - {3 + 2a)a - 3{2 + a)a^ + a^ f 



p,^k \ f3 J (1 + + a + 5aa - (5 - 3a)a^ - 3aa^ + Qa^) 

(78) 
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Figure 8: In the left panel we show s{a) = Cs/vb{a) and in the right panel Q^^^^/Q*^^, as 
functions of a, for the Jouguet detonation case. 
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Figure 9: Amplitude of the GW signal at the peak frequency from Jouguet detonations for 
(i/'H* = 10, 100, 1000, as a function of a (left) and as a function of vi,. The signal reaches 
a plateau at large a since it is bounded by the maximal possible value of the bubble wall 
velocity, Vb < I- 



which we plot in Fig. [91 

The peak amphtude depends quadratically on the duration of the phase transition, P 
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and grows steeply (also roughly quadratically) with a as long as a < 1. At a > 1 the 
dependence on a becomes rather weak, since the bubble velocity approaches its maximum, 
ff, — > 1. The maximum gravity wave density parameter which can be achieved by bubble 
collisions in a strongly first order phase transition a > 1, t'b — 1 and (5 ~ lOTi^, is of the order 
of 5 X 10"^'^. Note that our approximations break down for smaller values of /3, since we 
have neglected the factor ln(M/m) so that becomes the duration of the phase transition 
which must be significantly smaller than a Hubble time. 



5.3 GW from deflagrations 

In the case of defiagrations, Rif) = Vontt coincides with the shock front of the shock wave 
taking place in the symmetric phase. Therefore, we set Vout = "^^shock- The inner radius is the 
bubble wall, so that = Vh = V2- The fiuid velocity is given by 

vi = (79) 

1 - V^Vi 

where Vi is the incoming velocity of the symmetric phase fiuid into the front, in the rest 
frame of the discontinuity. 

In order to estimate the gravitational wave production, we follow the analysis of defia- 
grations presented in Appendix A of Ref. [20]. There, the authors numerically integrate 
the energy- momentum conservation equations at the defiagration front with different initial 
conditions. The velocity profile they find is in principle different both from the constant 
one occurring in planar deflagrations, and from the linear one that we took in our simpli- 
fled analysis. However, we choose two different cases which they have analyzed, for which 
the deflagration is quite strong and the velocity proflle is actually almost constant. In 
the flrst case, they set Vh = 0.1 and V{ = 0.09, and flnd f shock — 0.59 (and consequently 
s = "i^int/'^out = 'i^b/'i^shock = 0.17). In the second case, close to Jouguet deflagrations, they 
set Wb = 0.5 and V{ = 0.45, and flnd f shock — 0.73 (and consequently s = 0.68). The velocity 
proflles for these values are shown in Fig. 9 of [20] . Since in the case of deflagrations there is 
no simple relation among the velocity of the shock front and the parameter a denoting the 
strength of the phase transition, we simply leave Wshock as a free parameter in our analysis, 
for which we take the two values given above. 

The analysis of the shape of the spectrum is on the same footing as for detonations. The 
results of the unequal time correlator for the anisotropic stress are shown in Fig. [10] for flxed 
values of Ushock = 0.59 and Xc = 0.9n. The results for the incoherent and top-hat in the 
anisotropic stress cases are similar to the case of detonation. We recover the expected 
behaviour at large scales, and the amplitudes are comparable at low wave number. The peak 
is located at ^pcak — 4.2, corresponding to fcpcak — 4.2/5/fout = 4.2/i?(77fin) — 7r/i?(?7fin)- On 
the other hand, the coherent and top-hat in the velocity correlation cases are slightly different 
in the high frequency range, because the frequency of the oscillations is increased, due both 
to the smaller value of Vont = ^'shock = 0.59, and to the smaller value of s = 0.17. These 
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parameters in fact appear in the arguments of the Green functions in Eq. fl60l) and fl66!) . In 
the region where the power spectrum becomes negative, in the top hat in velocity case, the 
solution is no longer reliable. We therefore discard it and take as best approximation the 
unequal time correlator of the anisotropic stress power spectrum, as we did for detonations. 




Figure 10: Same as Fig.[6]but for the case of deflagrations. We fix Wghock = 0.59. Compared to 
detonations, the oscillations are enhanced in approximation Eq. (j60|) and (|66!) corresponding 
to top-hat in the velocity correlation. 



The result of Eqs. ( ITOHTTj) is very weakly dependent on the value of the velocity f shock- 
This is analogous to the weak dependence on a found for detonations. We do not display this 
dependence as it is very similar to the top left panel of Fig. [7| where the role of a is played by 
''^shock which is varied between 0.59 and 0.73. Correspondingly, when varying Xc, we obtain 
a dependence similar to the one shown in Fig. [71 For fixed = 0.9n and Wshock = 0.59 we 
can approximate the result as follows: 
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where int|dcfia(^) = ; — , 2^ , — rxs ' Z^ = 3.6. (80) 
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This is essentially identical to Eq. fl76l) except for the values of s and ff in ^^kin/^rad- 

case of deflagrations we cannot reduce the dependence exclusively to the two parameters a 
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Figure 11: The case of deflagrations: The left plot shows Q^^^^^/Q*^^ as a function of Wf. 
The red curve corresponds to the first case analyzed in [20], for which s = 0.17, with 
fb = 0.1 and f shock — 0.59 and the blue one to the second case with s = 0.68, = 0.5 and 
^^'shock — 0.73. The right plot shows the amplitude of the GW signal at the peak frequency 
given in Eq. (IHTj) as a function of V{ for different values of jS/TC,^. Solid lines are for s = 0.68 
corresponding to = 0.5 and f shock — 0.73, and dashed ones for s = 0.17 corresponding to 
t'b = 0.1 and I'shock — 0.59. As long as svf -C 1, we have f^kin/^rad °^ ^'f ^^'^ correspondingly 
dilP^^^h'^/dlnk oc vf. 



and (3, since no direct relation between the velocities of the shock front, the velocity of the 
bubble wall and the strength of the phase transition a is known in general. In this case, 
the velocity of the shock front is expected to depend also on the properties of the ambient 
plasma, like friction, and not only on the strength of the phase transition determined by a. 
The amplitude at the peak frequency is the same as the detonation formula fl77j) except that 
0.084 is replaced by 0.050. This factor accounts for the contribution of the integral at the 
peak frequency k ~ 4.2/3/fout (c/. Fig. [TOl) : 
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peak 



The first plot of Fig. [TT] shows ^kin/^rad ^ function of the fluid velocity ff for two values 
of s = fb/'i^shock, corresponding to the cases analyzed in [20]. Obviously, the value of Vf is in 
principle flxed once a, Wb and ?7shock are given; however, the relation among these parameters 
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is not known explicitly. Therefore, we have decided to keep Vf as a free parameter. The fluid 
velocity Vf induces the biggest variation in the amplitude of the signal, and is therefore the 
most relevant parameter determining ^kin/fi*ad- This appears clearly in the second plot of 
Fig. [m where we show the amplitude of the GW signal Eq. flSTl) as a function of ff , for the 
same two values of s and varying The dependence on s is negligible compared to the 

dependence on ff. 

6 Some comments on our approach 

In numerical simulations, the stochastic background of GWs arises from averaging over 
several deterministic realizations of bubble collisions. In our approach instead, to account 
for the intrinsic randomness of the nucleation process, we define the bubble velocity as a 
random variable. The source of the stochastic background of GWs is the tensor part of the 
anisotropic stress of the stochastic, homogeneous and isotropic velocity field. 

The calculation of the velocity correlation function in Sec. 13.2.11 implicitly assumes that 
the bubbles, and consequently the velocity configuration, are all independent from each 
other: the velocity correlation function is different from zero only if the points x and y are 
in the same bubble. Therefore, one may wonder in which sense this procedure is a model 
of bubble collisions. Indeed, in our model the only non-vanishing correlation coming from a 
(possibly) collided region is given by the sum of correlations from two independent bubbles: 
for every x and y we can find a bubble center position such that a bubble encompasses the 
two points, giving a non-zero result, provided that the probability to find a center in that 
position is not zero. This is accounted for by multiplying the correlation with the probability 
that a given point is in the broken phase at a given time, cf. Eq. fl28p . 

Therefore, even though we do not model in a deterministic way the resulting velocity field 
from the collision of bubbles, we do account for their overlap; it is the overlap from several 
bubbles which gives us a non-zero correlation function for the tensor part of the anisotropic 
stress, as opposed to the spherically symmetric situation described in Appendix A. The 
anisotropic stress correlation function, in fact, involves the four-point correlation function 
of the velocity field, cf. Eq. fll9l) . This quantity, in contrast with the two-point correlation 
function, is non-zero also for points x and y belonging to different bubbleqj. In other words, 
when we calculate the two-point correlation function of the velocity field Eq. fl28l) . we break 
spherical symmetry since every point in space becomes equivalent to another and there is 
no longer a center of symmetry. The center of symmetry xq has become a random variable, 
and we average over all its possible positions. Using the ergodic theorem, this is equivalent 
to an average over several realizations for the center positions, i.e. several realizations of the 
nucleation process, i. e. several realizations of the velocity field distribution. In this way, we 
account for the overlap of several bubbles, which breaks spherical symmetry and leads to a 
non-zero power spectrum of the tensor anisotropic stress, generating gravitational radiation. 

^{vaiyi)vbix)vc{y)vd{y)) D {va{x)vb{x)){vc{y)vdiy)). 
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7 Summary 



Since the GW signal from bubble collisions was already evaluated in Refs. [17-20], in this 
section we gather our final formulas and compare them with the ones given so far in the 
literature. We conform to the notation used previously, for the comparison to be straight- 
forward. The relevant formulas are compiled for instance in Ref. [24]: there, Eq. (4) gives 
the frequency at which the GW spectrum peaks and Eq. (5) shows the dependence of the 
amplitude of the GW signal at the peak frequency on the relevant parameters of the bubbles 
evolution. These formulas are valid in the case of Jouguet detonations, and do not show 
how the GW spectrum depends on other than the peak frequency. Before proceeding with 
the comparison, we recap the main assumptions in our calculation: 1) we assume radiation 
domination; 2) the fiuid velocity profile inside the bubble is linear; 3) Vi is our stochastic 
variable via the bubble center as we average over all possible positions of the bubble center. 
Therefore, we do not model collisions in a deterministic way. We rather account for all pos- 
sible configurations of the velocity field which include the case where bubbles overlap, even 
though we use the velocity profiles of uncollided bubbles; 4) we use the Wick theorem to 
express the 4-point correlation function in terms of the 2-point correlation function, even if 
the velocity is presumably not a gaussian random variable; 5) all bubbles have the same size 
R = v\y{ri — r^in); 6) we use the top hat approximation for the correlator of the anisotropic 
stress tensor at different times; 7) the enthalpy w and the Lorentz factor 7 do not depend 
on X. 

Using the notation 



,^n„„ (/) . "^1^ where / = ± and ^n,.,. - ""^^i^ 
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from the results given in Section [521 Eqs. ( 1761 1771) and Section [531 Eqs. (jHOl [8T1) . we find 
that h'^Q^aii{f) increases at small / as (compared to /^'^ in Kosowsky et al. [17]), and at 
large / it scales as f~^'^ (as in Kosowsky et al). The full spectrum is given by 
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Let us first discuss the peak frequency. It is given by 

/peak = -TT^ /Spoak - 4.5 where /5 = a*/3 and (85) 

Vont = Vh for detonations, Vont = ''^shock for deflagrations. (86) 

We remind that in our notations l3~^ expresses the duration of the phase transition: bub- 
bles are generated at the beginning of the phase transition and collide after a time given 
approximately by We have in fact set to 1 the logarithm relating /? and ?7fin — rji^, c.f. 
Eq. P6|) of section H] (this logarithm can be easily inserted in all the following equations: it 
will introduce a multiplying factor in the ratio (3/H and will also change the value of the 
integral (1601) ). Moreover, v^ut corresponds to either the characteristic bubble velocity (in 
detonations) or the shock wave velocity (in deflagrations). The characteristic frequency at 
the time of emission is /pcak/o* ~ P/vout- As expected, it is associated with the maximal 
size of the spherical fluid velocity configuration which generates the GW signal. This can be 
either the bubble itself, in the case of detonations, or the spherical shock front preceding the 
bubble, in the case of deflagrations. Since fshock < v^, the peak frequency for detonations is 
smaller than that for deflagrations (the size of the velocity configuration is bigger). Using 
= a^H^ and = a*/ao = To/T^(go/g^y^^, we can rewrite the peak frequency as 

/peak - — Which yields (87) 

271 Vont ri* U \g*J 



/p,,, ^ 1.12 X 10-2 mHz 



\im) 100 GeV Vont 



where of course /3/7i* = (3 / H^,. 

The above peak frequency is bigger than the one in Refs. [20,24] by a factor ~ 2/i;out- The 
factor l/fout5 that shifts the peak to higher frequencies for low velocities, is absent in [20,24], 
and is related to the fact that the characteristic frequency we find is determined by the size 
of the bubbles instead of the duration of the phase transition [4] . The property of causality 
of the source directly determines this characteristic frequency. Since the source is causal, the 
velocity correlation function goes to zero at the length-scale corresponding to the size of the 
bubbles. The same length-scale determines the peak in the velocity power spectrum, and 
from there it is transferred to the GW power spectrum [c.f. also the discussion in Ref. [43], 
where it is explained that gravity waves emitted during short cosmological events such as 
phase transitions typically inherit the wave number and not the frequency of the source) 0. 



^Note that, being in a cosmological context, we do not perform a time Fourier transform of the energy 
momentum tensor in contrast to Eq. (23) of Ref. [20]. In [20] the authors assume that "the frequency 
dependence of the spectrum is set by the timescale (as written before eq. (28) of [20]). Our finding, on 

the other hand, is that the frequency dependence of the spectrum is set by the bubble size R. Since in the 
detonation regime Vont approaches 1, in the simulations of Ref. [20] one does not really see this difference. 
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The factor 2, instead, comes from the details of our modehng of the source. It follows 
from the factor 4.5 in Eq. (l85l) (see Fig. [6] and Fig. [9] and discussion thereafter). We found 
that the velocity power spectrum has a characteristic wave number ~ 2.5/-R, corresponding 
quite well to that coming from the bubble diameter 2'k/2R. The time-dependent anisotropic 
stress power spectrum instead changes slope around 3/-R (see Fig. E]). The characteristic 
wave number of the GW spectrum is associated to -Rfin, the typical radius of bubbles at the 
end of the phase transition when bubbles collide, and is found to be /Cpeak ~ 4.5/-Rfln. This 
differs from the value obtained in Kamionkowski et al. /Cp^ak ~ 2/3 (see Fig. 7 of Ref. [20] and 
the associated uncertainty). For = 100, = 100 GeV and ~ 100, corresponding 

to a typical (first order) electroweak phase transition, we find /p^^k ~ 1 niHz/t>out to be 
compared with LISA's peak sensitivity which is estimated to be 2 mHz. The increase in 
the peak frequency that we obtain relative to Ref. [20] is actually welcome for probing the 
electroweak phase transition with LISAI^. 

Let us now discuss the peak amplitude. A simple order of magnitude estimate shows how 
the result depends on the duration and the energy density of the source. From the perturbed 
Einstein's equations 5G^u = ^txGT^u, one gets the following order of magnitude estimate for 
the amplitude of the tensor perturbation h (we drop indices for simplicity): 

(3^h^87rGT (89) 

where we inserted 1//3 as the characteristic time on which the perturbation is evolving, and 
T denotes the energy momentum tensor of the source. From Eq. f[T7|) and definition fISB]) . 
we can write 

T ~ Prad §^ . (90) 
^^rad 

We want to estimate the energy density in gravitational waves, defined in Eq. The above 
equation ( l89i) suggests that h ~ SnGT/P, and so we obtain 

Substituting Friedmann equation in the radiation dominated era = ^y^Prad, and con- 
sidering that the GW energy density evolves like radiation, we obtain for the GW energy 
density today the simple expression 

"-"-^"-©^(i:)^- '''' 

This shows that the GW energy density scales like the square of the ratio between the time 
duration of the source and the Hubble time, and the square of the energy density in the 
source. 



^This result should be taken with caution, see the "Note added" at the end of this Section. 
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More precisely, from Eqs. (177118 ip we get the following result: 



where 0.084 is replaced by 0.050 for deflagrations, and 



(94) 



with s = rint/-R: 



s = ^ , V{ = i^y^^ v\, = — — Y^:^ — — for Jouguet detonations, (95) 
s = "'^ for deflagrations. (96) 

There is no simple analytic relation between v^, V{ and Ushock in the case of deflagrations. 
For detonations, Vh{a) is given above, taken from Refs. [20, 34]. Using go = 3.36 and 
^rad/i^ = 4.2 X 10~^ we find the peak amplitude: 

/.^l],.., ^ 5.4 X 10- ^— ^ ( ^ ) (^) (^)\ (97) 
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where 9.8 is replaced by 5.8 for deflagrations. The signal vanishes if the fluid velocity is 
zero. It also vanishes when the shell's thickness is zero (s = 1), this happens if Vh = Cs 
(corresponding to Vf = 0) for detonations and for Vh = Wshock for deflagrations. Our result is 
to be compared with the expression reported in Eq. (5) of Ref. [24], which is only valid for 
Jouguet detonations: 

^^^^ - (l + a)2 + J [gj 

Here ^(a) = (0.715a + 4(3a/2)i/V27)/(l + 0.175a) is the parameter defined in Eq. (22) of 
Ref. [20], denoting the fraction of vacuum energy which goes into kinetic energy of the fluid 
(rather than thermal energy). Therefore, the combination a k is equivalent in our notation 
to the parameter ^^kin/^rad = Pkin/Prad- The above expression for n{a) is a fitting formula 
that the authors of [20] determined by integrating numerically the energy momentum tensor 
of the velocity profile corresponding to a Jouguet detonation, for different values of a. On 
the other hand, we have defined f^kin/^rad terms of the fluid velocity at the inner boundary 
of the velocity shell, starting from the approximated velocity profile given Eq. (!26|) . As it 
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should be, the dependence on a of the two parameters is comparable, see Fig. [12] where 
we used s = Cs/v^, V{ = (t>b — Cs)/(1 — fbCs) and Ub(Q!) = {cs + \/a'^ + 2a/3)/{l + a). In 
Fig. [ini we show the comparison between our peak amplitude fl98l) and the result used in [24] 
(Eq.EHD. 
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Figure 12: Comparison of an{a) defined in Ref. [20] (dashed line) with ^liri/^*a.d (solid line), 
where Vb{a) is given by Eq. (ff3l) . In the case of Jouguet detonations, both parameters are 
fully determined in terms of a. We remind that aK,{a) and ^kin/^rad defined exactly 

in the same way (see text), but they reflect the same physical quantity. 



flpeak h^/ [(H//?)2(100/g, ^^^^^ [(H/^)2(ioo/g. )V3] 




Figure 13: The case of Jouguet detonation. Comparison of the peak amplitude used in 
Ref. [24] (dashed line, Eq. (!99|) ) with our value (solid line, Eq. (1981) ) as a function of a in 
the left hand panel and as a function of Vb in the right hand panel. 

In summary, our amplitude for the signal is comparable (but can differ by one order 
of magnitude for some values of the parameters) to the one obtained with an inherently 
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Figure 14: The GW signal at the peak frequency, given in Eq. ( l98l) . as a function of the 
maximal value of the fluid velocity ff, for fixed s = 0.68. One can have Vf < Cg both for 
deflagrations and detonations, while Vf > Cs is possible only for detonations. One should 
keep in mind that in principle once s is fixed, Wf is not a free parameter. In the case of 
Jouguet detonations its value is well-known, Vi = Cs(l — s)/(s — c^); for deflagrations there is 
no analytic formula, and in order to derive Wf one needs to know as well as vi, the incoming 
fluid velocity in the frame of the bubble discontinuity (c./. section [3l2|) . Nonetheless, the 
aim of this figure is just to show the order of magnitude of the GW signal when the fluid 
velocity is taken as a free parameter. This gives a reasonable estimate, since the dependence 
on s is small, c.f. Fig. [TT] 



different method, that of numerical simulations in the envelope approximation [17-20]. This 
confirms that the details of the collision's modeling are not so crucial and what really matters 
at the end is the size of the velocities involved in the process. Note also that we computed 
the GW energy density spectrum without using the Weinberg formula i.e. without making 
the wave zone approximation. The latter assumes that the observer is at a distance much 
larger than the dimension of the source, while our source is spread over the entire universe. 

Kamionkowski et al. obtain a non-vanishing signal in the limit of vanishing thickness. 
This can be understood as follows. In Ref. [17-19], they first study bubble collisions taking 
place in vacuum, so that the source of GW was the energy momentum tensor of the scalar 
field: being given by the spatial gradient of the scalar field, this is non-zero only at the bubble 
wall. In their later work Ref. [20] , they consider bubble collisions in a thermal bath and in this 
case they use as the GW source the energy-momentum tensor of the relativistic fluid rather 
than that of the scalar field. Nevertheless, they keep using the envelope approximation. This 
is why they obtain a large signal even if the shell of fluid velocity has vanishing thickness. 
In contrast, by construction, there is no signal in this limit in our model, meaning that the 
kinetic energy of the fluid (which is our only source as we did not include the gradient energy 
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of the scalar field) is different from zero only over a finite volume. On the other hand, we 
can extrapolate our results to the deflagration regime. Typically, for the detonation case, in 
the limit fb Cg Ref. [20] finds a non-zero GW signal, since in the envelope approximation 
this is just the lower bound of the bubble wall velocity. Within our model, in this limit, the 
thickness of the non-zero velocity shell goes to zero and the GW signal as well; however, 
this only shows the break down of the detonation regime, and the necessity of treating the 
problem in the deflagration approach. 

In Figure [TU we plot the peak amplitude as a function of the maximal fluid velocity V{. 
In contrast with the case of Jouguet detonations, for deflagrations, we cannot express the 
signal as a function of a only. Indeed, for small velocities, the relation between a and the 
bubble velocity will depend on the interactions between the bubble wall (the Higgs field) and 
the particles in the thermal bath. In any given model of a first-order phase transition, one 
can in principle compute the bubble wall velocity and the consequent fluid velocity profile, 
see e.g. Ref. [42] for the case of a weakly first-order EW phase transition, and Ref. [44] for 
the strongly first-order phase transition presented in Ref. [6]. Provided that the released 
latent heat is large, one can obtain a large bubble wall velocity. In addition, one would 
have to look carefully at the physics of deflagrations to derive the relation between the fluid 
velocity and that of the bubble wall, see for instance Ref. [35]. 

To conclude, our main new results can be summarized as follows: 

• Our description applies to both detonation and deflagration (the fluid velocity is non- 
zero over a finite volume rather than on an infinitely thin wall and v ^ 1). 

• Our peak frequency is parametrically larger (oc 1/v). 

• We provide an analytic expression for the shape of the spectrum. 

The possibility that the signals discussed here, if they are produced at the electroweak 
phase transition, could be detected with LISA will be discussed in an upcoming publica- 
tion [29]. We confirm that the GW signal coming just from bubble collisions is observ- 
able only for very large fluid velocities. Indeed, EISA's best sensitivity is not much below 
Qh^ ~ 10~^^. Such a value can be reached if ~ 10 and V{ ~ 0.2. More realistic values 

j3/H^, ^ 10^ would require V{ ~ 0.5 in order to lead to an observable signal. As mentioned 
in the introduction and as will be presented in [29], there are other sources (magnetic field 
and turbulence) of GW during phase transitions and the signal from bubble collisions is just 
one contribution. 

Note added 

Our results for the position of the peak and the slope of the spectrum at large frequencies 
strongly depend on the time structure of the anisotropic stress. In our modelization, it grows 
with time and is suddenly switched off at the end of the transition. This discontinuity is 
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somehow unphysical. Indeed, our Eq. [29] does not take into account the fact that towards 
the end of the transition, there is no more intersecting region but only one single bubble. 
As discussed in details in [45], if the source is switched off smoothly and we consider the 
coherent case as indicated by numerical simulations [46] and if the relevant correlation length 
is taken to be the size of the uncollided region rather than the bubble size, then the peak 
position does not depend on the velocity and is given by kpeak ~ P- Besides, numerical 
simulations [46] which are carried out under the thin-wall approximation also indicate that 
the slope of the spectrum at high frequencies typically scales as 1/k. 



Acknowledgments 

We thank Riccardo Sturani for reading the manuscript. CC and RD acknowledge support 
by the Swiss National Science Foundation. 



A One single bubble 

If there was only one bubble, there would not be any gravitational radiation since a single 
bubble corresponds to a spherically symmetric distribution of energy and momentum which 
cannot emit gravitational waves. In this Appendix we show that the anisotropic stress Hij 
from a single bubble is purely scalar: this means that there exists a function / such that 



SijA)f , 



(100) 



where Tjj denotes the energy momentum tensor of a bubble. Such a scalar component is 
projected out by the tensor projection operator which is given in Fourier space in Eq. (jS]), 
and does therefore not contribute. In this case, from the above equation the function / must 
satisfy the condition 

AV = ^d,d0^ . (101) 

We now demonstrate that this condition is always satisfied for a single, spherically symmetric 
bubble. For one spherically symmetric bubble we have, up to the constant enthalpy p + p, 



Vi{r)vj{r) - -6ijv'^{r) 



so that 



didj 



Vi{r)vj{r) - -SijV^ir) 



vv 



vv" + v'^ + 5 h 



2r2 



4 , 2 2 
-V V H — V 
3 r 



(102) 
(103) 

(104) 

(105) 
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For the third equal sign we have used spherical symmetry, we work with the Ansatz v = 
v{r)er, and the prime denotes the derivative w.r.t. r. On the other hand, for a spherically 
symmetric function / we have 



Therefore, if we find a function / which satisfies 



6 r 6 \ r J 



or, equivalently 




(106) 



the anisotropic stress (11021) is equal to the expression given in Eq. fllOOl) . therefore, it is a 
pure scalar. This is indeed the case, since Eq. (11061) is an ordinary linear differential equation 
which always has a solution for a given radial velocity v. 

Hence a single spherically symmetric bubble only generates scalar perturbations and 
does not contribute to the tensorial part of the energy momentum tensor Iljj, which is the 
source of gravity waves. This is of course not surprising given that spherically symmetric 
configurations only have scalar degrees of freedom. Nevertheless, we have added this brief 
calculation here since one might draw wrong conclusions from the fact that the anisotropic 
stress of a single, spherically symmetric bubble does not vanish. 

Another way to arrive at the same result is to show that the tensor projection operator 
given in Eq. ([8]), 



vanishes when applied to a spherically symmetric stress tensor. For a single bubble one has 
then 



B Calculation of the velocity field power spectrum 

In this appendix we explain how to evaluate the velocity correlation function Eq. ( 128|) and 
Fourier transform it to obtain the velocity power spectrum Eq. fl38p . The quantity we need 
to calculate is the tensor 



MijW = PikPjl — -;j^PijPki 



with 



(107) 
(108) 



liij = MijkiTki = Mikji {vk{r)vi{r)) = . 




(109) 
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The intersection volume Vi varies with the distance between x and y, as shown in Fig. [3l 
with r = |x — y|. We choose an orthonormal basis with 62 || x — y. One identifies four 
different regions: setting a=|x — y|/2 = r/2, they are given by the hmiting values 



< a < 



R - rint 



rint <a< 



2 ' 

R + ^'int 



R - nnt 

2 

R + ^int 



<a< Tint , 

< a< R. 



;iio) 



The intersection volume is symmetric under rotations around 62, so in order to perform the 
integral in fll09p we choose cylindrical coordinates with z \\ 62 and p \\ 63. We have 

d^Xo = pdpdz dip , xq = (pcos(v9), 2;, psin(</))) , x=(0,— a,0), y = (0, a, 0) . (Ill) 



Substituting the above formulas in (11091) and performing the integration in dip, one sees that, 
because of cylindrical symmetry, the tensor ly is diagonal. Evaluating (11091) reduces simply 
to calculate the two integrals 



hi = /33 = TT / dzdpp 

'Vi 



l22 = 2n ! dz dp p {z'^ - a^) . (112) 
Jvi 



The limits of integration, here generically denoted with V^, depend in fact on the variable a. 
As an example, in the region < a < (i? — rint)/2, the first integral becomes: 



h 



21 rpi 

IT I dz dp p^ 

Z2 Jo 

Zl = r-^at + a , Z2 = R — a. 



ppi 
dz I dp p^ 



Zl rpi rZ2 rpi 

dz dp p^ + dz dp p^ 
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Pi = VR' - (« - zr 

P4 = - (« + zy 



P2 



\lrtt-{a + zf , Pa = ^Jrl^,-{z-ay 



Analogous expressions can be found for the remaining three regions, and similarly for /22. 
The final result of the integrations are the two continuous functions In (a. Tint, -R) and 
122(0, Tint, -R), as a function of the variable < a < i?, which are too long expressions 
to be written explicitly here. 

Knowing In = I^^ and /22, we impose the condition of statistical homogeneity and 
isotropy for the velocity field, meaning that we impose the tensorial structure 



lij = f Sij + g f-if-j = In + {!'. 



22 



hi) riTj 



(113) 



where the second equality is a straightforward consequence of our choice f = x — y || 62. 
The two-point correlation function for the velocity field takes the final form (c/. Eq. [321) 



,2 r 



(t;,(x,t)t;,(y,t)) = 0(t)-^ 



h 



1 A I 22 



- I 



11 



(114) 
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where the pre-factor is independent of a = r/2 and the volume Vc is given in Eq. fl30|) . 

In order to know the power spectrum, we have to Fourier transform the above equation 
with respect to the variable r. Because of homogeneity and isotropy, we get a delta function 
in momentum. We rewrite the Fourier transform of the second term in the sum, which is 
direction dependent, in terms of derivatives with respect to the wave vector, and we obtain: 



Vr 



1 22 — III 



;ii5) 



The Fourier transform of a function only of r gives a function only of wave number k. Know- 
ing this, we can re-express the partial derivatives, to obtain expression flM|l and foUowings. 
The Fourier integrals, as for example 



47r 



drr sin (/cr 



Iii{r, r;^t,R) 
Vrir) 



(116) 



need to be further divided in the sum over the four integrals corresponding to the regions 
described above, depending on the value of r, since in each of these regions the integrand 
takes a different form. Again, we do not write the complicated full expression of the result, 
for which we found the fitting formulas fl39l) and fHOj) shown in Fig. |H 



C Large and small scale part of the GW power spec- 
trum 

We have seen that the gravitational wave power spectrum, independently of the different time 
approximations, always grows like at scales larger than the peak scale k~^^^ 2^ i?(?7fin)/4.5. 
The reason for this general behaviour is that the source is uncorrelated at these scales: 
therefore, the power spectrum of the anisotropic stress source is simply the incoherent sum 
of uncorrelated regions, and is white noise. The white noise behaviour for the anisotropic 
stress in turn determines the k^ increase for the GW power spectrum. 

On the other hand, we also saw that for the small scale part of the spectrum we always 
recover roughly a decrease, with the only exception of the completely incoherent case 
(c/. Eqs. ( !62][64|) ). The small scale decrease can be understood from dimensional arguments. 
In this Appendix we present general arguments for the origin of the above mentioned power 
laws for the large and small scale part of the GW power spectrum. 

Let us first concentrate on the large scale limit. We start with a generic velocity power 
spectrum showing a peak at a characteristic scale k = L~^. The scale L, corresponding in 
our case to the bubble diameter L = 2R, may depend on time. We assume that the large 
and small scale behaviors are given by two power laws. 
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satisfying the conditions n > —3, m < —3 so that the energy density is dominated by the 
contribution at the peak. The pre-factor denotes the average energy per unit enthalpy 
of the source. For bubbles, we have n = for the function A{k) and n = 2 for B{k), cf. 
Eqs. (inniinD and the discussion thereafter (we will treat the small scale decrease of these 
functions later on). The anisotropic stress power spectrum is given by the convolution of the 
velocity power spectrum (c/. Eq. [24l) . Setting Q = qL and K = kL, and neglecting angular 
dependencies, we have 

/»oo r /»1 /"oo 

Il{k)(x dqq^P{\k-fi\)P{q)(xv^L' / dQQ^^+^PHK - Q\) + dQQ""-^^ P{\K - Q\) 
Jo Uo Jl 

(118) 

In the large scale limit k <^ L ^, in the second integral we can safely neglect K with respect 
to Q and simply set P(|K — Q|) ^ Q"^. In the first integral we have to be a little more 
careful. If n < —3/2, the main contribution to the integral comes from the divergence at 
Q — s> K and the integral picks a behaviour i^2n+3 ^^j, physical bubble situation, we have 
n = 0, n = 2; we therefore choose to analyze only the case n > —3/2, for which we find: 



U{k ^ 0) oc v^L^ 



dQQ 



2n+2 



dQQ 



2m+2 



1 



2n 



2m + 3 



;ii9) 



Therefore, the anisotropic stress at large scales is white noise. Going back to Eqs. ([6]) and 
( fT5l) . we see that for a constant anisotropic stress power spectrum, the GW power spectrum 
behaves like 



dn{k) 



d\n{k) 



oc k^\h'{k)\^ oc eU{k) 



(120) 



Here we have re-expressed the conformal time derivative appearing in Eq. ([6]) in terms of 
a derivative with respect to x = kr], and we have taken into account that converting the 
double integral over the Green function in Eq. f|T5l) into a double integral with respect to 
conformal time, induces an additional factor k"^. We therefore recover the observed large 
scale behavior oc k^. 

We now turn to the small scale limit, k ^ L~^. In this case, our fits to the functions 
A{k) and B{k) decrease like k~'^. As argued in Section [3l3l this power law is transferred to 
the small scale behaviour of the anisotropic stress power spectrum, cf. Eqs. f H2]H3|) . which 
takes the form 



(n,,-(k, r)n*,(q, r)) cx 5(k - q)i?(r)3X(Alr)) 



:i2ii 



where T{K) decreases like i^"^ for K ^ 1, and we neglect all the other time dependent 
factors, which are irrelevant for the argument presented here. Inserting the above expression 
in Eq. (fT5|) . one finds 



\h'{k)\'^^^ jdyj dzR' 



V) cos(?/ 



z)X{K{t)) oc ^ 



dy / dzy\os{y-z)X{y) , (122) 
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where for the second equahty we have used the fact that R{t) oc r. Substituting the above 
formula into Eq. ([6]) rewritten in terms of the derivative with respect to x, one finally obtains 



dfl{k) 



d\n{k) 



oc oc dy j dz y'' cos {y - z)I{y) , (123) 



where the double integral simply causes a modulation in the spectrum. We recover, in the 
large wave number limit, the /c~^ decrease, which is a simple consequence of dimensional 
analysis since X is a function of y = krj only in the relevant regime. This power law usually 
acquires small corrections due to the fact that the integral in Eq. (11231) is not completely 
independent of k. Depending on the approximation for the unequal time anisotropic stress 
power spectrum, we actually found power laws /c~^ with j3 in the range 1.8 < (3 < 2.2. Within 
our approximation, these may well also be logarithmic corrections to the slope /5 = 2. 

In the case of the totally incoherent unequal time approximation, the argument is changed 
because of the delta function in time: the form of the anisotropic stress power spectrum is 
now (c/Eq. [MD 

(n,,-(k, r)n;.(q, 0) oc 6{k - q)5(r - C)i?(r)2j(i^(r)) . (124) 

Repeating the above argument, one easily finds that the large wave number power law is 
now changed to k~^ {k~^'^ with the corrections from the integral). 
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